Następnie, za pomocą testów zgodności, zostanie dopasowany najlepszy rozkład do częstości strat oraz do ich dotkliwości. Zostanie wygenerowanych 20000 scenariuszy i zostanie wyznaczony symulacyjny rozkład strat. Procedura wyznaxczania OpVaR i OpES zostanie powtórzona 200 razy. Został mi przydzielony plik22.csv.
library(dplyr)
library("tidyverse")
library(fitdistrplus)
library(actuar)
library(DT)
library(gridExtra)
library(BatchGetSymbols)
options(scipen = 0, digits = 4)
dane <- read.csv("plik22.csv")
czestosc <- dane %>%
group_by(linia) %>%
count(rok)
wartosci <- unique(dane$linia)
ramka_danych <- data.frame("Year"=unique(dane$rok))
for(i in 1:length(wartosci)){
x <- czestosc[czestosc$linia==wartosci[i],3]
ramka_danych <- cbind(ramka_danych,x)
}
colnames(ramka_danych) <- c("Year",wartosci)
datatable(ramka_danych)
dwu1 <- fitdist(ramka_danych[,2], distr="binom",
start = list(prob=0.5),
fix.arg = list(size=13))
dwu2 <- fitdist(ramka_danych[,4], distr="binom",
start = list(prob=0.5),
fix.arg = list(size=15))
dwu3 <- fitdist(ramka_danych[,7], distr="binom",
start = list(prob=0.5),
fix.arg = list(size=27))
uj1 <- fitdist(ramka_danych[,3], method="mme",distr="nbinom")
uj2 <- fitdist(ramka_danych[,5], method="mme",distr="nbinom")
uj3 <- fitdist(ramka_danych[,6], method="mme",distr="nbinom")
uj4 <- fitdist(ramka_danych[,8], method="mme",distr="nbinom")
ramkapvalue <-
data.frame(c(gofstat(dwu1)$chisqpvalue,
gofstat(dwu2)$chisqpvalue,
gofstat(dwu3)$chisqpvalue,
gofstat(uj1)$chisqpvalue,
gofstat(uj2)$chisqpvalue,
gofstat(uj3)$chisqpvalue,
gofstat(uj4)$chisqpvalue))
colnames(ramkapvalue) <- c("P.value")
rownames(ramkapvalue) <- c("Buss_Distr", "Damage", "External_Fr",
"Com_ban", "Empl_pract", "Execut_Del",
"Internal_Fr")
datatable(ramkapvalue)
############################################################
Buss_Distr <- dane[dane$linia==wartosci[1],3]
Com_Ban <- dane[dane$linia==wartosci[2],3]
Damage <- dane[dane$linia==wartosci[3],3]
Empl_Pract <- dane[dane$linia==wartosci[4],3]
Execut_Del <- dane[dane$linia==wartosci[5],3]
External_Fr <- dane[dane$linia==wartosci[6],3]
Internal_Fr <- dane[dane$linia==wartosci[7],3]
lista <- list(Buss_Distr,
Com_Ban,
Damage,
Empl_Pract,
Execut_Del,
External_Fr,
Internal_Fr)
pvalue <- list()
for(i in 1:length(lista)){
r1 <- fitdist(lista[[i]], "weibull")
r2 <- fitdist(lista[[i]], method="mme", "lnorm")
r3 <- fitdist(lista[[i]], method="mme", "exp")
r4 <- fitdist(lista[[i]], method="mme", "gamma")
r5 <- fitdist(lista[[i]], distr ="pareto", start = list(shape = 1, scale = 500))
test <- gofstat(list(r1,r2,r3,r4,r5),
fitnames = c("weibull", "lnorm", "exp", "gamma","pareto"))
chisqpvalue <- ifelse(test$chisqpvalue < 0.05, "rejected", "not rejected")
pvalue[[wartosci[i]]] <- data.frame("KSTest"=test$kstest,
"ADTest"=test$adtest,
"ChisqTest"=chisqpvalue,
"ChisqPvalue"=test$chisqpvalue)
}
## Warning in cov2cor(varcovar): 'diag(.)' miał 0 lub wpisy NA; nieskończony
## rezultat jest wątpliwy
## Warning in sqrt(diag(varcovar)): wyprodukowano wartości NaN
## Warning in sqrt(1/diag(V)): wyprodukowano wartości NaN
## Warning in cov2cor(varcovar): 'diag(.)' miał 0 lub wpisy NA; nieskończony
## rezultat jest wątpliwy
datatable(t(as.data.frame(pvalue)), options = list(pageLength = 4))
param_buss <- fitdist(lista[[1]], "weibull")
param_com <- fitdist(lista[[2]], method="mme", "gamma")
param_dmg <- fitdist(lista[[3]], "weibull")
param_empl <- fitdist(lista[[4]], "weibull")
param_exe <- fitdist(lista[[5]], "weibull")
param_ext <- fitdist(lista[[6]], "weibull")
param_int <- fitdist(lista[[7]], method="mme", "lnorm")
wyk1 <- ggplot(data = data.frame(lista[[1]]), aes(x = lista[[1]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_buss$estimate[1],
param_buss$estimate[2])) +
labs(title = "Buss_distr", x = "Buss_distr") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk2 <- ggplot(data = data.frame(lista[[2]]), aes(x = lista[[2]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Gamma"),
fun = dgamma,
args = list(param_com$estimate[1],
param_com$estimate[2])) +
labs(title = "Com_ban", x = "Com_ban") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Gamma")) +
theme_minimal()
wyk3 <- ggplot(data = data.frame(lista[[3]]), aes(x = lista[[3]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_dmg$estimate[1],
param_dmg$estimate[2])) +
labs(title = "Damage", x = "Damage") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk4 <- ggplot(data = data.frame(lista[[4]]), aes(x = lista[[4]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_empl$estimate[1],
param_empl$estimate[2])) +
labs(title = "Empl_Pract", x= "Empl_Pract") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk5 <- ggplot(data = data.frame(lista[[5]]), aes(x = lista[[5]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_exe$estimate[1],
param_exe$estimate[2])) +
labs(title = "Execut_Del", x="Execut_Del") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()+
xlim(c(0, 200))
wyk6 <- ggplot(data = data.frame(lista[[6]]), aes(x = lista[[6]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_ext$estimate[1],
param_ext$estimate[2])) +
labs(title = "Extern_fr", x="Extern_fr") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk7 <- ggplot(data = data.frame(lista[[7]]), aes(x = lista[[7]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Log-norm"),
fun = dlnorm,
args = list(param_int$estimate[1],
param_int$estimate[2])) +
labs(title = "Intern_fr", x="Intern_fr") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Log-norm")) +
theme_minimal()
liczba_zdarzen <- 20000
liczba_powtorzen <- 200
VaR <- c()
ES <- c()
for(i in 1:liczba_powtorzen){
Buss_czest <- rbinom(liczba_zdarzen, size = 13, prob = dwu1$estimate)
Dmg_czest <- rbinom(liczba_zdarzen, size = 15, prob = dwu2$estimate)
Ext_czest <- rbinom(liczba_zdarzen, size = 27, prob = dwu3$estimate)
Com_czest <- rnbinom(liczba_zdarzen,
size=round(uj1$estimate[1],0),
mu=uj1$estimate[2])
Empl_czest <- rnbinom(liczba_zdarzen,
size=round(uj2$estimate[1],0),
mu=uj2$estimate[2])
Exec_czest <- rnbinom(liczba_zdarzen,
size=round(uj3$estimate[1],0),
mu=uj3$estimate[2])
Int_czest <- rnbinom(liczba_zdarzen,
size=round(uj4$estimate[1],0),
mu=uj4$estimate[2])
ramka_czestosc <- data.frame("Buss"=Buss_czest,
"Com"=Com_czest,
"Damage"=Dmg_czest,
"Empl."=Empl_czest,
"Exec."=Exec_czest,
"Ext."=Ext_czest,
"Int."=Int_czest)
straty_ramka <- data.frame()
for(i in 1:nrow(ramka_czestosc)){
strata_buss <- sum(rweibull(ramka_czestosc[i,1],
shape=param_buss$estimate[1],
scale=param_buss$estimate[2]))
strata_com <- sum(rgamma(ramka_czestosc[i,2],
shape=param_buss$estimate[1],
scale=param_buss$estimate[2]))
strata_dmg <- sum(rweibull(ramka_czestosc[i,3],
shape=param_dmg$estimate[1],
scale=param_dmg$estimate[2]))
strata_empl <- sum(rweibull(ramka_czestosc[i,4],
shape=param_empl$estimate[1],
scale=param_empl$estimate[2]))
strata_exe <- sum(rweibull(ramka_czestosc[i,5],
shape=param_exe$estimate[1],
scale=param_exe$estimate[2]))
strata_ext <- sum(rweibull(ramka_czestosc[i,6],
shape=param_ext$estimate[1],
scale=param_ext$estimate[2]))
strata_int <- sum(rlnorm(ramka_czestosc[i,7],
meanlog = param_int$estimate[1],
sdlog = param_int$estimate[2]))
straty_ramka <- rbind(straty_ramka,
c(strata_buss,
strata_com,
strata_dmg,
strata_empl,
strata_exe,
strata_ext,
strata_int)
)
}
colnames(straty_ramka) <- colnames(ramka_czestosc)
straty_ramka$Suma <- rowSums(straty_ramka)
wartosci <- sort(straty_ramka$Suma, decreasing = T)
kwantyl <- 0.999
VaR <- c(VaR, quantile(wartosci, kwantyl))
indeks_do_wycięcia <- length(wartosci) * (1-kwantyl)
wyc_wartosci <- wartosci[1:indeks_do_wycięcia]
ES <- c(ES, mean(wyc_wartosci))
}
datatable(ramka_czestosc)
## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html
datatable(straty_ramka,
options = list(
columnDefs = list(list(
targets = "_all",
render = JS("function(data, type, row, meta) {",
"if(type === 'display' && typeof data === 'number') {",
"return parseFloat(data).toFixed(2);",
"} else {",
"return data;",
"}",
"}")
))
))
## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html
datatable(
data.frame("VaR"= c(mean(VaR), sd(VaR), sd(VaR)/mean(VaR)), "ES" = c(mean(ES),sd(ES), 100*sd(ES)/mean(ES)),
row.names = c("Średnia","Odchylenie st.", "Współczynnik zmienności (%)"))
)
\[ E_t = A_tN(d_1) - Le^{-r(T-t)}N(d_2), \]
\[ d_1 = \frac{\ln(A_t/L) + (r + (\sigma^2)/2)(T-t)}{\sigma\sqrt{T}} \]
\[ d_2 = d_1 - \sigma\sqrt{T} \]
\(N(-d_2) = 1 - N(d_2)\).
\[ \left\{ \begin{align*} E_0 &= V_0N(d_1) - De^{-rT}N(d_2), \\ \sigma_{E}E_0 &= N(d_1)\sigma_{V}V_0 \end{align*} \right. \]
first_date <- as.Date("2023-01-01")
last_date <- as.Date("2023-12-31")
l.out <- BatchGetSymbols(tickers = c("ASOZY"),
first.date = first_date,
last.date = last_date,
do.cache = FALSE)
ceny <- l.out$df.tickers$price.open
zwrot <- diff(log(ceny))
liczba_dni <- length(zwrot)
zmiennK <- sd(zwrot)*sqrt(liczba_dni)
wartKap<-83000303*mean(ceny)
r<-0.0625
T<-1
zobow<-9765000000
#Funkcja licząca błędy pomiędzy wartością kapitału oszacowaną oraz rzeczywistą wraz ze zmiennością kapitału
MertonModel<-function(x){
V0<-x[1]
zmiennV<-x[2]
d1<-(log(V0/zobow)+(r+zmiennV^2/2)*T)/(zmiennV*sqrt(T))
d2<-d1-zmiennV*sqrt(T)
A<-V0*pnorm(d1)-zobow*exp(-r*T)*pnorm(d2)-wartKap
B<-pnorm(d1)*zmiennV*V0-zmiennK*wartKap
return(A^2+B^2)
}
# Funkcja minimalizująca błędy
min_blad <- function(x){
wynik <- optim(c(V0=x[1],zmiennV=x[2]),MertonModel)
return(c(wynik$par[1], wynik$par[2]))
}
# Prawdopodobieństwo bankructwa
PD <- function(V, zob, r, zmienV, T){
d1<-(log(V/zob)+(r+zmienV^2/2)*T)/(zmienV*sqrt(T))
d2<-d1-zmienV*sqrt(T)
return(pnorm(-d2))
}
wyniki <- min_blad(c(11000000000,0.0))
V0<-wyniki[1]
zmiennV<- wyniki[2]
pd <- PD(V0, zobow, r, zmiennV, T)
datatable(as.data.frame(data.frame(data=c(V0, zmiennV, pd), row.names = c("Wartość aktywów", "Zmienność", "PD"))))